% input for sample discrete ordinates code

%clear
format short

p = 4;
if ( p==0 )
    xcm    = [ 0  2  4  6  8  10 ];        
    xfm    = [  1  1  1  1  1    ]*10;
    src    = [  1  1  1  1  1    ];   
    mt     = [  4  4  4  4  4   ];
elseif ( p==1 )
    % sample 1, its = 91, 34, 48
    xcm    = [ 0  2  4  6  8  10 ];        
    xfm    = [  1  1  1  1  1    ]*100;
    src    = [  1  1  1   1  1   ];   
    mt     = [  4  2  4  3  4    ];
elseif (p==2)
    xcm    = (0:20)/200;  
    xfm    = ones(1,20)*10;
    src    = [ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
    mt     = [ 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4];
elseif (p==3) %
    xcm    = [0 6 12 18 24 30];
    xfm    = [ 6 6  6  6 6];
    src    = [1 1  1 1 1];
    mt     = [5  5  5 5 5];
elseif (p==4) % homogeneous, highly-scattering slab                                
%       SI  |  CMR  | CMFD  | CMSC  | mesh(mfp)
%      910  |   62  |diverge|  69   | 5
%      910  |   27  |diverge|  35   | 2.5   
%      910  |   23  |    29 |  34   | 1.25   (stable for up to ~1.67mfp)
%      910  |   56  |    11 |  92   | 0.625   
%      910  |diverge|     7 | 344   | 0.3125    
%
% Note, for 0.625 mfp mesh and c=0.9999, CMFD still takes just 7 its for the
% 1e-6 convergence on phi whereas unaccelerated run takes 14441 its.
%
%     xcm    = 0:5:100;                   
%     xfm    = 32*ones(1,length(xcm)-1);
%     xcm    = 0:2.5:100;
%     xfm    = 16*ones(1,length(xcm)-1);    
%     xcm    = 0:1.25:100;                  
%     xfm    = 8*ones(1,length(xcm)-1);
%     xcm    = 0:0.625:100;                 
%     xfm    = 4*ones(1,length(xcm)-1);
    xcm    = 0:0.3125:100;                 
    xfm    = 2*ones(1,length(xcm)-1);
    src    = ones(1,length(xcm)-1);
    mt     = 4*ones(1,length(xcm)-1);
end



data   = [ 1.00   0.50   0.50    % mid scatter
           1.50   1.50   0.00    % pure absorber
           0.50   0.30   0.20    
           1.00   0.01   0.99    % high scatter
           1.00   0.00   1.0  ]; % pure scatter   

for a = [0]
input   =   struct(       ...
    'numg',            1, ...     % number of groups
    'numm',            1, ...     % number of materials
    'xcm',           xcm, ...     % slab bounds
    'xfm',           xfm, ...     % number of fine meshes
    'mt',             mt, ...     % slab material ids
    'data',         data, ...     % mat comp's
    'src',           src, ...     % volume source
    'ord',            32, ...     % number of ordinates
    'maxit',        1e6, ...     % max iterations
    'maxerr',       1e-6, ...     % max pointwise phi error
    'adj',             0, ...     % adjoint flag
    'bcL',             0, ...     % left boundary condition
    'bcR',             0, ...     % right boundary condition
    'acc',             a  ...     % acceleration (0=none,1=cmr,2=cmfd
    );
hold on

ntime = 1;
time = zeros(ntime,1);
for t = 1:ntime
    tic
    if (input.acc==0)
        [phi,psi,x,J,er0]   = sn_one_d(input); c = 'k';
    elseif (input.acc==1)
        [phi,psi,x,J,er1]   = sn_one_d(input); c = 'b--';
    elseif (input.acc==2)
        [phi,psi,x,J,er2]   = sn_one_d(input); c = 'c-.';
    else
        [phi,psi,x,J,er3]   = sn_one_d(input); c = 'g:';
    end
    time(t) = toc;
end
disp([' average over 10 runs is ', num2str(mean(time)), ' +/- ', num2str(std(time))])

% a sampling of phi from various cm's
[ phi(1) phi( 2 ) phi( 3) ]

end

plot(x,phi,c,'LineWidth',2),grid on
% xlabel('x [cm]'), ylabel('\phi(x) [n/cm^2]'), title('Sample 1-D, 1-G Problem')
% axis([0 10 min(J)*1.1 max(phi)*1.1])
% grid on
% %phi=[mean(phi(1:10)) mean(phi(11:20)) mean(phi(21:end))]
% %plot(psi)
% 
% figure(2)
% semilogy(1:length(er0),er0,'k',1:length(er1),er1,'r',...
%         1:length(er2),er2,'b','LineWidth',2)
% legend('none','cmr','cmfd')
% grid on

